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SUMMARY 

The  response  of  an  aircraft  as  a function  of  time  on  encountering  an 
isolated  ramp  gust  is  derived  from  its  response  to  a unit  step  gust.  Two 
FORTRAN  programs  are  described  which  treat  separately  straight  ramp  and  smooth 
ramp  (one-minus-cosine)  gust  profiles.  The  (critical)  gust  length  causing 
maximum  d}rnamic  response  is  determined  and  responses  to  simple  gust  patterns 
(pairs)  are  investigated. 
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I INTRODUCTION 

The  problem  discussed  in  this  Report  is  the  following: 

Given  the  response  F(t)  of  an  aircraft  flying  at  a certain  speed  to  a 
unit  step  gust,  predict  the  response  to  a family  of  ramp  gusts  and  determine  the 
critical  gust  length  at  which  the  extreme  response  occurs. 

The  investigation  was  treated,  purely  as  a mathematical  and  computational 
exercise;  the  soundness  and  applicability  of  the  physical  theory  {eg  Jones S was 
taken  for  granted.  Jones  recalls  in  that  Report  that  the  earliest  airworthiness 
requirements  were  based  on  response  to  discrete  gusts  but  that  recently  more 
emphasis  has  been  laid  on  irregular  turbulence,  the  implication  being  that 
responses  can  be  deduced  adequately  from  stationary  random  process  theory.  He 
then  argues  that  large  aircraft  loads  cannot  be  satisfactorily  predicted  by 
spectral  analysis  and  that  the  consideration  of  discrete  ramp  gusts  is  preferable 
if  extreme  responses  are  to  be  reliably  estimated.  Effectively,  the  argument  is 
that  the  power  spectrum  (PSD)  method  estimates  the  root  mean  square  response 
whereas  one  is  really  interested  in  the  extreme  behaviour  and  that  the  non- 
Gaussian  character  of  turbulence  can  be  so  marked  as  to  undermine  the  usefulness 
of  the  rms  value.  In  particular,  the  PSD  method  is  not  to  be  trusted  when  the 
dominant  aircraft  mode  is  well  damped. 


The  problem  is  the  familiar  one  that  the  weight  of  interest  falls  in  the 
tail  of  a distribution  where  probability  levels  are  only  poorly  defined  by 
empirical  evidence.  Any  extrapolation  from  the  relatively  well  defined  body  of 
the  distribution  is  unsafe  unless  independent  information  on  the  shape  of  the 
distribution  is  available.  In  addition,  Jones  points  out  that  dynamic  response 
is  produced  by  ohangee  in  wind  velocity  so  that  one  must  determine  the  probability 
distribution  of  velocity  gradients.  Many  power  spectrum  investigations  have 
only  considered  the  velocity  distribution. 


To  cope  with  situations  in  which  non-Gaussian  velocity  gradients  are 

significant,  Jones  has  developed  a statistical  discrete  gust  model  wherein  the 

intensity,  w of  a ramp  gust  is  related  to  its  gradient  distance  (H)  by 
H 


w 


H 


(1) 


where  w^  is  derived  directly  from  measurements  and  varies  with  altitude.  This 
dependence  of  gust  intensity  on  gradient  length  is  valid  over  the  range  of  wave- 
lengths (X)  for  which  the  Kolmogoroff  (proportional  to  X^  ; X 2H)  power 
spectrum  of  atmospheric  turbulence  applies.  Fortunately,  the  Kolmogoroff  spectrum 

f 
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appears  to  be  valid  far  beyond  the  wavelengths  associated  with  extreme  aircraft 
response. 

The  ramp  (see  Fig  1)  may  be  either  straight  or  'one-minus-cosine'  (smooth). 

Wj(x)  - 0 f X < H (2) 

W2(x)  - - cos  ; 0 $ X « H . (3) 

dw  dw 

The  latter  has  the  advantage  that  — « 0 at  x “ OtH  (whereas  — is 

dx  dx 

discontinuous  at  the  end  points)  and  is  generally  preferred. 

In  the  next  section  we  derive  the  ramp  response  from  the  unit  step  response 
and  the  gust  profile.  Later  sections  explain  the  computational  details  and  the 
determination  of  that  gust  length,  H which  causes  the  worst  response,  y • 

This  is  extended  in  section  7 to  deduce  the  most  adverse  gust  pair. 

Although  the  problem  appears  straightforward  at  first  sight,  being  basically 
one  of  numerical  integration  and  interpolation  there  are  complicating  features. 

In  addition,  although  the  physical  model  is  only  approximate,  we  did  not  wish  to 
further  degrade  it  with  rough  calculations.  Our  efforts  to  ensure  efficiency, 
accuracy  and  stability  in  the  numerical  operations  are  duly  described. 

Thus,  in  section  3 we  describe  the  integration  procedure  based  on  cubic 
spline  fitting  of  the  step  response,  explaining  the  considerable  economies  which 
are  possible  with  this  formulation.  The  straight  ramp  generates  a discontinuity 
in  the  slope  of  the  ramp  response  which  can  be  most  troublesome  (if,  for  example, 
the  'corner'  is  smoothed)  if  not  properly  handled  as  described  in  section  4. 

Extra  trial  gust  lengths  are  needed  to  refine  the  estimate  of  the  critical  gust, 

H , and  more  especially,  the  critical  time,  t . These  are  kept  to  a minimum  by 
an  efficient  interpolation  procedure  (see  section  5). 

The  end  result,  we  believe,  are  computer  programs  which  can  determine  the 
critical  ramp  gust  and  gust  pair  with  errors  which  are  negligible  compared  to 
chose  arising  from  Che  fundamental  physical  assumptions  and  idealisations. 

2 DERIVING  THE  RAMP  RESPONSE  FROM  THE  STEP  RESPONSE 

2. I Assumptions 

The  important  assertion  which  motivated  Che  present  work  is  that  the  res- 
ponse to  an  isolated  ramp  gust  of  a given  length  and  profile  shape  is  a meaning- 
ful quantity  to  calculate.  Jones ^ claims  that  such  a representation  is  usually 
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reasonable  at  least  in  the  vicinity  of  the  critical  gust  length  (i<?  that  isolated 
gust  which  generates  the  worst  behaviour).  In  other  words,  despite  turbulence 
being  comprised  of  gusts  of  all  lengths,  only  a limited  range  is  important  and 
the  magnitude  of  the  worst  response  can  be  fairly  well  predicted  by  identifying 
the  critical  isolated  gust  (or,  more  generally,  sequence  of  gusts). 

A second  major  assumption  is  that  the  system  is  linear.  A corollary  is 
then  that 

Y(Wq,H)  - WqYO.H) 

where  y is  the  (extreme)  response. 

The  linearity  assumption  also  facilitates  consideration  of  gust  pairs 
whereby  overswing  tuning  can  result  in  a greater  response  than  any  single  gust 
would  cause.  The  effect,  of  course,  is  unimportant  for  well  damped  modes.  Gust 
pairing  is  therefore  an  extension  of  the  basic  concept  to  moderately  damped  modes 
and  involves  finding  the  critical  gust  length,  , which  results  in  the  worst 

overswing  as  well  as  that  (H)  which  gives  the  extreme  primary  response.  The 
most  adverse  situation  is  when  the  two  gusts  occur  with  opposite  sign  and  appro- 
priate spacing  (see  Fig  2 and  section  7)  so  that  the  first  overswing  due  to 

the  first  gust  coincides  with  the  primary  peak  due  to  the  second  gust.  The 
spatial  extent  of  the  critical  gust  pattern  is  therefore 

H-  - ii  + H-  + H . (4) 

2 Os 

More  complex  sequences  of  gusts  are  not  considered.  However,  it  is  possible  to 

extend  the  concept  to  include  very  lightly  damped  modes  and  a recent  paper  by 
2 

Jones  envisages  the  superposition  of  responses  to  as  many  as  eight  successive 
gusts.  Of  course,  the  probability  of  the  critical  pattern  occurring  decreases 
as  the  number  of  component  gusts  is  increased. 

2.2  Theoretical  considerations 

Let  the  step  response  be  F(t)  and  the  ramp  response  be  ((i(H,t)  . 
Essentially,  the  calculation  of  involves  the  numerical  integration  of  F(t)  . 
Note  that  strictly  speaking  we  should  write  F(v,t)  but  it  is  customary  to 
regard  the  aircraft  velocity  as  constant  during  an  encounter  with  turbulence. 

It  is  not  a simple  (to  linear)  matter  to  determine  (Ji  (and  hence  y(H))  at 
different  speeds  - F(t)  must  be  recalculated. 
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A ramp  gust  can  be  thought  of  as  a series  of  n step  gusts  (see 
Fig  3).  Thus  we  have  to  combine  together  up  to  n step  responses,  each 
scaled  and  appropriately  time  shifted.  The  linearity  assumption  is  essential 
if  the  means  of  combination  is  to  be  simple  addition.  The  response  at  time 
t*  will  be  given  by 


j ■ min  (n,t'/6t)  . Letting  6t  -►  0 


where  n6t 


Or,  making  the  transformation  t ••  t 


F(t)dt 


■ 0 when  t < 0 or  t > H/v  . Expressed 
equations  (2)  and  (3)  become 


The  limits  reflect  the  assumption  that 
as  functions  of  time  and  setting  w.  > 1 


Hence  (7)  gives,  respectively 


F(t)dt 


max(0,t'-H/v) 


(t’  - t)F(t)dt 


max(0,t 


7 


Thus  the  straight  ramp  response  at  time  t'  is  merely  the  area  under 
the  step  response  between  t'  and  some  earlier  time,  either  zero  or  t'  - H/v  . 
Furthermore,  dit ferent iat ing  CM  gives 


3H 


«^,(H,t') 


0 


(II) 


(12) 


From  equation  (II)  we  can  envisage  three  possible  situations  concerning 
the  time,  t^  , at  which  the  maximimi  response,  y(H)  occurs.  As  depicted  in 
Fig  4 can  be  greater  than,  less  than,  or  equal  to  H/v  . The  first  two 

present  no  computational  difficulties  because  * 0 but  the  third  case 

could  be  troublesome  because  the  discontinuity  in  coincides  with  t^  . 

Now,  in  any  numerical  scheme  integrals  can  only  be  computed  at  discrete 
(usually  pre-determined)  points.  If  one  seeks  maxima  and  minima  of  a function 
defined  only  as  a set  of  points  then  sixne  interpolation  is  necessary.  Of 
necessity  the  interval  between  adjacent  points  will  be  represented  by  a fully 
continuous  interpolating  fvinction  and  uifficulties  can  be  expected  if  the 
original  function  is  less  well  behaved. 

Trouble  was  indeed  encountered  with  the  straight  ramp  when  the  derivative 
discontinuity  at  t ■ H/v  fell  within  the  same  integration  interval  as  t 

Y 

The  resulting  interpolation  errors  superimposed  a ripple  on  the  y(H)  curve 
making  impossible  the  accurate  determination  of  H . The  problem  only  arises  if 
F(0)  0 but  unfortunately  this  is  comnx^n  because  F(t)  is  the  response  to 

a step  gust.  The  counter-measures  adopted  are  described  in  section  4.  There 
was  no  corresponding  difficulty  with  the  one-minus-cosine  ramp  because 
9i>,j/3t  is  everywhere  smooth. 

It  can  also  be  surmised  that  t will  increase  with  H until  the  first 
zero,  ty  , of  F is  reached.  For  H > vt^^  one  would  expect  t_^  “ t^^ 

(t»>’  case  2 ot  Fig  41  and  >110  to  decrease  as  H . F.quat  ion  U2'>  provides  the 
condition  to  be  satisfied  at  H . First  we  observe  that  t ^ H/v  is  impossible. 
For  t > H/v  we  have 

Y - I H*F(7  - H/v)  . 
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But  from  ( I t ) > 


F(t)  - F(t  - H/v) 


3 -V  - 
i H^F(t) 


If  t - H/v  then  neither  derivative  (11)  or  (12)  is  separately  zero  but  we 
still  have 


giving  y 


j H*F(H/v)  . 


Hence  equation  (14)  remains  valid  when  t * H/v  . It  is  useful  for  checking 
purposes  since  the  computer  program  does  not  use  it  directly  but  instead  locates 
H by  interpolation. 

Equation  (9)  for  the  straight  ramp  gust  encounter  can  be  solved  very 
economically  because  the  integrand  is  independent  of  H . One  merely  has  to 
store  once  and  for  all  the  accumulated  area  under  F(t)  at  predetermined  time 
points  and  then  use  appropriate  bandwidths  when  considering  different  H . 
Equation  (10)  is  far  more  complicated  because  the  step  response  is  weighted  by 
an  H-dependent  factor.  Furthermore,  practical  checking  formulae  analagous  to 
(13)  and  (14)  cannot  be  given. 

3 EVALUATION  OF  INTEGRALS 

The  first  operation  in  both  the  straight  ramp  and  the  smooth  ramp  programs 
is  to  fit  a cubic  spline  to  F(t)  . The  two  spare  degrees  of  freedom  are  taken 
up  by  third  derivative  continuity  at  the  second  and  penultimate  points  (knots). 
Subroutine  TB04A  of  the  Harwell  Subroutine  Library  is  employed  for  this  purpose. 
A particular  virtue  of  the  cubic  spline  is  that  the  area  between  knots  is  easily 
calculated  and  will  be  accurate  to  third  order  even  when  the  abscissae  are 
unequally  spaced  (af  Simpson's  rule). 


- ''Ui> 


where  F'  ■ dF/dt  . The  integral  in  (9)  is  therefore  readily  evaluated  while 
that  in  (10)  requires  further  manipulation  as  fellows. 


r 


Haintl  tho  cubic  spliuc,  tbc  'wciubCcd*  (in  the  sense  ot  (10))  area  between 
twc<  adjacent  knots  is  given  by 


/ 


sin  a(t'  - t)K(t)dt 


- r5^oos 


L 


cos  n ( t ' - t ) 


2 

sin 

vt(t' 

- 1) 

It* 

. iLiLa 

3 

cos 

a(t  ' 

- t) 

A 

K"’(t) 

4 

sin 

a(l  ’ 

- t) 

(lb) 


If  we  now  consider  the  total  integral  over  ni  intervals,  the  fact  that 

K.K’  and  K"  are  continuous  at  the  knots  and  that  K" ' is  constant  on  each 

interval  causes  nK*st  terras  to  cancel  leaving 


. r , T ™ 

I sin  o(t'  - t)K(t)dt  “ *'  cos  o (t ' - t ) . . . ♦ ~ y (sina(t'~  *i-i^ 


0 i-l 

- sin  o ( t ' - t ) ) . 


(W) 


Now,  from  equation  (10)  we  have  a ■ nv/H  , t “ t*  . Hence  sin  a(t'  ~ ’ 


cos  a(t'  - i ) - I while  at  the  lower  integration  limit  w«'  must  distinguish  the 
m 


two  cases  (i)  ~ • (it)  ^0  * ' 


In  the  first  case  we  have: 


(i)  sin  o(t'  - t^^)  - 0 , cos  a(t'  “ *^))  * ~ ' • 


Therefore  the  first  term  on  the  right  hand  side  of  (17)  beconH>s 


]L(i'-)_t.>I(lJ  - H/v)  _ F"(t*)  y*'(t’  - H/v) 


nv/H 


(KS) 


(nv/H)' 


Regarding  the  second  term  on  the  right  of  (17),  all  intervals  except 
possibly  the  first  (depending  whether  or  not  H/v  is  an  integral  multiple  of 
At)  will  be  of  the  satin'  sire.  At  . I'herefore, 
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m ID 

I-I-r- 

i-l  i-2 


(sin  a(t'  - t._|)  - sin  a(t'  - t^))  + (sin  a(t'  - t^) 


sin  a(t ' - 1 1)) 


where  Fl"  is  the  third  derivative  of  F(t)  in  the  ith  interval  following  t^  . 
But  t^  “ t'  - H/v  , tj  - t^  - Atj  (say)  , t'  - tj^  “ (m  - i)At  and 

m ■ — rounded  up  if  necessary. 
vAt 


Therefore  ~ 

i“l  i“2 

m 

- 2 sin  X cos  ^ (m  - i + J)At)-  F”' 


in  aAt , 


cos  — (m  - X + 
H 


J)At)-  F”' 


H ■ 


In  the  second  case  we  have: 

(ii)  All  intervals  are  of  length  At  and  the  integrand  does  not  vanish  at  rhe 
lower  limit.  The  first  term  on  the  right  of  (17)  becomes 


F(t')  F"(t')  F(0)  . F"(0)  F'(0)  . 

— ^ ^ ^ cos  at  + r; — cos  at ^ — sin  at 

a 3 a 3 2 

a a a 


while,  following  the  above  analysis,  the  sum  reduces  to 


m 

2 sin  ^ Fl"  cos  (m  - i + })At 


where  m ■ t* /At  . 

It  is  therefore  necessary  to  evaluate  F"  at  each  knot  and  F”'  on  each 
interval  and  to  store  them.  In  addition,  in  order  to  evaluate  efficiently  the 
integrals  for  all  t'  at  given  H it  is  helpful  to  first  store 


cos  (j  - i) 


for  j ■ 1,2  ...  integer  part  of 


The  only  quantities  which  then  have  to  be  recalculated  at  each  time  point  are 
cos  at'  and  sin  at'  (for  use  in  (20)).  When  t'  > H/v  even  this  becomes 
unnecessary  although  F(t'  - H/v)  and  F"(t'  - H/v)  may  then  have  to  be 
specially  interpolated  from  the  cubic  approximation  over  the  relevant  interval 
for  use  in  ( 18) . 

Potential  trouble  spots  in  the  above  formulation  are  the  heavy  dependence  on 
third  derivatives  and  the  presence  of  F' (0)  and  F"(0)  in  (20)  because  no 
restrictions  are  placed  on  the  derivatives  generated  at  the  end  points  of  the 
spline.  Tests  with  analytic  step  responses  have  verified  that  the  numerical 
scheme  can  be  very  accurate  but  ultimately  of  course,  in  any  practical  case  much 
will  depend  on  the  nature  of  F(t)  and  the  sampling  interval,  At  (see  also 
section  8) . Perhaps  the  more  important  feature  of  the  comparison  with  analytic 
expressions  was  the  validation  of  the  programming  of  equations  (18)  to  (21). 

4 OBTAINING  y(H) 

The  previous  two  sections  explained  how  the  ramp  response  is  calculated  at 
equally  spaced  time  points  for  each  trial  gust  length,  H . We  now  have  to 
detect  its  primary  peak  and,  where  applicable,  the  maximum  overswing.  No 
assumptions  can  be  made  about  the  shape  of  the  curve  which  makes  the  location  of 
global  extrema  all  the  more  difficult.  The  basic  procedure  in  both  programs  is 
to  search  the  set  of  ii>(t)  for  the  largest  positive  and  negative  values  and  to 
then  interpolate  y(H)  and  t^  using  d(|i/dt  . The  situation  depicted  in 
Fig  5 will  therefore  be  wrongly  analysed  as  shown  but  it  is  hoped  that  At  can 
be  specified  sufficiently  small  to  render  such  a shortcoming  unimportant. 

As  noted  in  section  2,  when  using  the  straight  ramp  model  we  have  to  be 
particularly  careful  to  evaluate  ())(H/v)  explicitly  because  d<|i/dt  is  there 
discontinuous.  At  other  absicissae  d(j>/dt  can  be  found  directly  from  the  tabu- 
lated step  response  using  (11).  We  therefore  have  to  consider  the  following 
possibilities. 

(1)  4)(H/v)  is  the  largest  ordinate  found  initially, 

then  either  (a)  d<p/dt  changes  sign  at  H/v  in  which  case  t^  = H/v  , or 

(b)  d((i/dt  does  not  change  sign  at  H/v  , te  t^  lies  between  H/v 

and  one  or  other  of  the  adjacent  abscissae. 

(2)  <|>(H/v)  is  not  the  largest  ordinate  found  initially  but  jAt  is.  Then  t^ 

lies  somewhere  between  (j  - l)At  and  (j  + l)At  . Again,  if  H/v  lies  within 

this  interval,  care  must  be  taken  to  locate  the  nearest  point  to  jAt  which 
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brackets  the  peak  ao  that  t and  y(H) 
t*>  t^  lies  in  one  of  the  intervals  ((j 
(jAt,H/v). 


can  be  accurately  interpolated, 

- OAt.iAt),  IjAt.li  ♦ DAt),  (H/v.iAt), 


The  process  of  interpolating  the  peak  is  rot.<rely  that  of  finding  a stationary 
point  of  a cubic  given  the  function  value  and  first  derivative  at  two  bracketing 

abscissae. 


The  smooth  ramp  model  differs  in  that  t - H/v  does  not  present  any  special 
difficulties  but  on  the  other  hand,  obtaining  d<i>/dt  is  not  so  straightforward. 
In  fact  differentiating  (10)  gives 


2 2 
s V 

2ir 


/ 


niax(0,t'-H/v) 


cos  ^ - OF(t)dt  . 


{22) 


instead  of  evaluating  the  above  integral  we  have  elected  to  fit  a cubic  spline 
through  interpolate  the  required  extrema  from  it  having  located 

the  maximum  and  minimum  tabulated  points  by  inspection  as  in  the  straight  rami' 
mode  1 . 

5 FINDING  ii,  y(H)  AND  t^(H) 

Both  the  primary  peak  of  «Kt)  and  the  largest  overswing  are  determined  for 
each  trial  value  of  H supplied  by  the  user.  Actually,  we  seek  the  extreme 
positive  and  negative  responses  y^(H) , Y_(H)  . It  was  originally  suggested  that 
we  call  the  larger  (in  absolute  value),  and  the  smaller,  y_  but  tltis 

seemed  to  be  nx're  confusing  and  could  (if  • ****  have 

defined  them,  intersect)  complicate  the  shape  of  Y(10  causing  discontinuities 
In  dy/dH  and  corresponding  difficulty  in  locating  II  (see  Fig  h) . Using  the 
revised  convention,  >(H)  is  nx're  likely  to  be  unimodal  and  only  after  and 

Y have  been  estimated  do  we  decide  which  is  the  primary  peak  and  which  is  the 
overawing.  It  is  not  necessarily  the  case  that  the  sign  of  the  primary  is  the 
same  as  that  of  the  first  peak  in  the  ramp  response  nor  that  the  worst  ('verswing 
occurs  after  the  primary. 

In  order  to  find  H it  is  essential  that  the  user  svipplies  trial  values  of 
H which  bracket  H , or  nx're  precisely,  that  the  extreme  occurs 

neither  at  the  smallest  nor  the  largest  gust  length  supplied.  In  priiu'iple  the 
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program  could  select  its  ovm  trial  values  of  H , choosing  gust  lengths  anywhere 

between  vAt  and  vt  (where  t is  the  last  point  of  the  given  step 

max  max 

response).  In  practice  it  was  felt  that  the  response  to  specific  gusts  could  be 

of  interest  so  the  programs  expect  a minimum  of  three  trial  gust  lengths.  The 

user  should  avoid  H < vAt  or  H > vt  and  should  be  careful  to  supply  enough 

max 

time  history,  F(t)  . Strictly  speaking,  sufficient  step  response  should  be  input 
to  enable  y(2H)  to  be  calculated  but  this  may  not  always  be  possible  (see 
section  6) . 


Assuming  these  conditions  have  been  fulfilled  we  will  get  a best  estimate  of 
H and  lower  and  upper  bounds  and  . We  then  use  a safeguarded  quadratic 

interpolation  scheme  outlined  below  to  successively  improve  H until 
max(H  - H,  H - H, ) is  sufficiently  small.  Several  further  ramp  responses  will 
have  to  be  calculated  but  we  feel  that  the  labour  is  justified.  It  may  be 
thought  possible  to  economise  calculation  by  only  evaluating  f in  the  vicinity 
of  t^  but  unfortunately  t^  can  change  discontinuous ly  (and  hence  t cannot 
safely  be  interpolated)  as  the  following  example  shows. 


Consider  the  step  response  shown  in  Fig  7a  in  conjunction  with  the  straight 

Then  (see 


ramp  model  and  let  area  A3  exceed  area  Al  and  t^^  - tQ2  ^ 


Fig  7b)  t ■ H/v  until 
Y 


H 


vt 


increasing  to  t^^  when 


01 
1 

v(t 


is  reached  after  which 


t ~ 

Y 01 


until  H is 


sufficiently  large  that  A3b  »-  Al  when  it  switches  di  scont  inuous  ly  to  t ^ 


03  " *^02^ 


Eventually  t will  revert  to 


01 


at 


where 


14 

*^14 

/ 


F(t)dt  - Al 


14 


H/v 


This  is  a powerful  argument  in  favour  of  computing  extra  ramp  responses  rather 
than  attempting  to  interpolate  t directly  from  the  responses  at  the  given  values 
of  H . Although  Y and  H may  perhaps  be  satisfactorily  derived  from  the 
initial  set  it  is  apparent  that  t could  be  grossly  in  error  and  hence  the  most 
adverse  gust  pair  predicted  (see  section  7)  would  be  totally  misleading. 


We  now  describe  the  quadratic  approximation  procedure  used  to  improve  H . 

The  simplest  algorithm  would  be  one  which  always  maintained  three  points  bracket- 
ing the  peak  so  that  interpolation  was  always  possible.  However,  if  the  curve 
was  very  asymnetric  such  a process  would  be  inefficient.  For  example,  in  Fig  8, 
point  3 would  always  be  retained  and  slow  convergence  would  be  experienced  because 
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of  its  remoteness  from  H . A preferable  scheme  is  to  maintain  wherever  possible 
the  three  highest  points  even  if  it  means  extrapolating  H . An  artificial  bound 
would  then  be  needed  to  safeguard  the  extrapolation.  Such  algorithms  are  common- 
place in  one-dimensional  search  routines  appearing  in  optimisation  programs  and 

3 . 4 . 
are  discussed  by  Brent  and  Gill  and  Murray  . If  H is  the  highest  point  so  far 


found  and  and 

is  defined  by 


are  the  current  bounds  on  H then  the  artificial  bound 


H + 0(H^  - H) 


H + B(Hy  - H) 


H ^ p 


H < p 


where  p - J + H^)  ; B - 1(3  “ . 

If  the  extrapolation  predicts  a point  outside  then  is  used  instead. 

This  strategy  becomes  comparable  with  the  usual  bracketing  technique  as  the 
skewness  of  the  curve  decreases.  In  fact,  the  programs  work  with  the  curve 
Y(log  H) , making  the  determination  of  H even  easier,  t is  not  interpolated 
but  is  obtained  directly  from  (|>(H,t)  as  of  course  is  y . 

A further  refinement  prevents  the  evaluation  of  at  values  of  H which  are 

too  close  together.  The  user  supplies  a parameter  specifying  the  relative 
accuracy,  AH/H  , required  and  the  program  locates  H such  that 

max(log^(H/Hj^)  , log^(Hy/H))  < AH/H  (24) 


while  keeping  the  minimum  separation  between  any  pair  of  (logarithmic)  trial 
AH 

gust  lengths  > . 

6 SENS IT  IV m FACTOR 

Jones*  defines  A , the  gust  length  sensitivity,  a measure  of  the  sharpness 
of  the  peak  in  y(H)  by 

1 /2y(H)  - y(2H)  - Y(H/2)>f 


log^,2 


2itY(R) 


An  alternative  expression  which  ignores  gusts  longer  than  H is 


w \ tiy(h)  / 

These  quantities  are  evaluated  for  both  Y^(*0  'r_(H)  • 


( 
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Note  that  the  calculation  of  y(2H)  may  require  the  provision  of  an  uneconomic 
amount  of  step  response.  Likewise,  the  user  should  beware  of  suggesting  too  short 
a trial  gust  length  or  of  occasions  where  H is  so  short  that  H < 2vAt  . The 
programs  will  ignore  trial  gust  lengths  < vAt  but  will  attempt  to  determine 
y(H/2)  in  similar  circumstances  so  that  X may  be  evaluated.  Note  that  for 
ramp  gusts  of  very  short  gradient  distance  (straight  or  smooth)  we  have  the 
asymptotic  result 

(^(H.t)  ~ H^F(t)  as  H ->  0 . (27) 

Thus  a unit  step  gust  can  effectively  be  regarded  as  being  equivalent  to  a ramp 
gust  of  length  1 ft.  However,  the  programs  do  not  exploit  this  feature, 
ie  the  maximum  and  minimum  of  F(t)  are  not  determined. 

7 GUST  PAIRS 

Having  calculated  identify  the 

primary  peak  with  the  larger  of  y^  ^nd  |y_|  , the  smaller  becoming  the  worst 
overswing.  We  Chen  deduce  the  most  adverse  gust  pair  situation  as  follows. 

The  maximum  possible  response  is  (y^  + 1y_1)  occurring  at  t^  or  t_ 
whichever  is  the  later.  Let  us  assume  for  clarity  that  t_  > t^  . Then,  as 
noted  in  section  2,  the  worst  gust  pair  is  a ramp  gust  of  length  H_  followed 
by  one  of  opposite  sign  of  length  H^  , separated  by  H^  where  (see  Fig  2) 

H = v(t  - t ) - ii  . (28) 

s - + - 

It  is  conceivable  that  H^  < 0 will  result  but  again  we  leave  aside  the  inter- 
pretation of  such  a situation. 

8 ACCURACY 

As  stated  at  the  beginning  of  this  paper,  we  have  tackled  a specific  mathemati- 
cal problem  leaving  it  to  others  to  assess  Che  results  and  the  validity  of  the 
models  used.  The  interaction  with  our  computer  programs  really  only  occurs 
through  the  specification  of  AH/H  . The  fixing  of  this  parameter  depends 
primarily  on  one's  assessment  of  the  accuracy  of  the  discrete  ramp  model.  Having 
made  such  an  appraisal  it  is  then  up  to  the  user  to  ensure  that  the  step  response 
is  sufficiently  well  represented  for  the  desired  accuracy  to  be  attained  for  it 
must  be  remembered  that  the  numerical  scheme  has  a finite  accuracy  and  spline 
fitting,  integrating  and  interpolating  are  all  sources  of  error.  The  question 
then  is  - how  close  is  the  problem  wliich  the  program  has  accurately  solved  to 
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that  which  one  has  tried  to  specify?  In  a practical  situation  the  only  insight 
available  is  by  rerunning  the  program  with  fewer  points  defining  the  step  response 
and  comparing  results.  If  they  are  too  different  then  a finer  tabulation  of 
P(t)  is  probably  necessary. 

Both  the  straight  ramp  and  smooth  ramp  programs  have  been  tested  on  a two 
parameter  family  of  analytic  step  responses  suggested  by  Jones  and  which  may  be 


written  as 


where  ■ I - 


F(t)  “ e ^^fcos  nt  + 


a » I ; 0 < e < > 


sin  Sit 


With  At  ••  0.5  and  v - 100  , H and  t were,  with  the  smooth  ramp  program 
always  obtained  to  within  5Z  while  y was  even  better  (<1Z).  However,  the 
dominant  factor  affecting  numerical  accuracy  is  the  density  of  tabulated  points 
relative  to  the  scale  on  which  F varies.  The  period  of  the  above  function 
exceeds  2it  by  definition  of  S)  . At  “ 0.5  is  thus  not  a large  spacing  by  any 
means.  Tests  have  suggested  that  the  straight  ramp  program  is  relatively  more 
accurate  probably  because  of  the  need  to  spline  fit  in  the  smooth  ramp 

program. 

Normally,  F(t)  would  be  obtained  by  solving  a set  of  differential  equations 
in  which  case  At  would  automatically  reflect  the  variability  of  F so  that 
hopefully  no  appreciable  loss  of  detail  need  occur  in  feeding  F(t)  to  the  ramp 
gust  programs.  A possible  snag  is  that  differential  equation  solvers  tend  to 
vary  the  step  size  thereby  tabulating  F(t)  in  a manner  unacceptable  to  the 
smooth  ramp  program  which  requires  equally  spaced  points. 

In  sunmary,  we  are  confident  that  <1Z  error  in  H,  t and  y is  attainable 
on  an  ICL  1900  series  computer  for  arbitrary  F(t)  , providing  F(t)  is  well 
tabulated  (ta  we  expect  the  numerical  methods  to  be  able  to  handle  many  more 
than  the  200  data  points  currently  allowed  by  the  dimensions  of  arrays). 

Obviously,  precision  will  suffer  when  the  programs  are  run  on  a machine  possessing 
lower  floating  point  accuracy  and  double  precision  working  may  then  become 
compulsory. 

9 USING  THE  PROGRAMS 

There  should  be  little  difficulty  implementing  the  programs  because  they 
have  been  written  in  standard  FORTRAN  and  the  data  required  is  minimal.  The 
following  quantities  need  to  be  input  on  channel  3. 
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(a)  TMAX  the  time  of  the  last  point  in  the  step  response.  Set  negative  if 

analytic  F(t)  (programmed  in  FUNCTION  FUNC(T))  is  to  be  used. 

(b)  TINT  the  time  interval  at  which  F(t)  is  to  be  calculated  or  supplied  on 

channel  1.  In  the  smooth  ramp  program,  TINT  is  necessarily  also  the 
interval  at  which  the  ramp  response  is  evaluated.  If  F(t)  is  supplied 
in  tabular  form  to  the  straight  ramp  program  then  TINT  is  not  read  and 
the  time  points  are  expected  on  channel  1 and  can  be  arbitrarily 
spaced.  This  facility  will  probably  be  little  used  (but  see  section  8) 
because  the  smooth  ramp  model  is  preferred  and  the  same  step  response 
data  cannot  be  used  inanediately. 

(c)  VELY  the  speed  of  the  aircraft  in  ft/s. 

(d)  DT  the  time  interval  at  which  the  straight  ramp  response  is  to  be  cal- 

culated. Not  read  in  the  smooth  ramp  program  because  DT  - TINT 
necessarily. 

(e)  NH  the  number  of  trial  values  of  H supplied  (NH  ^ 3) 

(f)  H(I),  I - l,NH  trial  gust  lengths.  H to  be  bracketed; 

t > H/v  > At  preferably  (see  section  5). 

IDdX 

(g)  HLTOL  accuracy  required,  AH/H  . 

(h)  IPR  print  control  parameter  giving  different  levels  of  output  on  channel  2 

as  follows 


IPR 


and 


for  trial  H(I)  and  H 


worst 


for  all  extra  gust  lengths  invest i- 


- 1 data  summary;  y 
gust  pair: 

- 0 as  IPR  - - 1 plus  > and 

gated  including  JH  and  2H  . Also  the  integrals  in  (9)  or (10)  are 
output  at  intervals  of  At  for  H ■ H and  JH  . 

■ I as  IPR  ■ 0 plus  the  tabulated  spline  fit  to  the  step  response. 

> 2 as  IPR  " 0 plus  ramp  responses  for  each  H(I)  supplied  by  the  user. 

The  multiplying  factor  vH  or  vH  is  also  output  for  con- 

verting Che  integrals  to  4i(H,c)  . 

■ 3 as  IPR  ••  1 and  IPR  2 combined. 

(i)  lANAL  Set  to  zero  if  subroutine  ANALYTIC  to  evaluate  (l>(H,t)  analytically  is 
either  not  provided  or  not  being  used.  Otherwise,  input  any  non-zero 
integer  value. 

This  completes  the  channel  3 data.  If  TMAX  > 0 then  the  step  response  must 
be  tabulated  on  channel  1.  Both  programs  currently  include  coding  for  F(t)  given 
by  (29)  and  the  straight  ramp  version  can  determine  iti(H,t)  by  integrating 


analytically.  They  expect  a and  5 to  be  input  immediately  after  TMAX  (<  0) . 
Figs  9 to  12  typify  the  computer  output  (IPR  » 3).  The  straight  ramp  program  pro 
duces  similar  output  except  that  the  columns  of  the  spline  fit  headed  D2F/DT2  and 
D3F/DT3  are  absent  because  they  are  not  needed. 

10  CONCLUSICWS 

Ac  the  request  of  J.G.  Jones  of  Flight  Systems  Department,  RAE,  Bedford, 
two  FORTRAN  programs  have  been  written  to  facilitate  the  identification  of  Che 
discrete  gusts  which  produce  the  greatest  aircraft  response.  The  input  to  the 
programs  is  the  calculated  response  as  a function  of  time  to  a unit  step  gust. 
From  this  one  can  determine  Che  behaviour  on  encountering  a gust  whose  intensity 
profile  is  of  ramp  form.  The  process  is  in  principle  quite  straightforward  but 
converting  the  theory  into  an  algorithm  suitable  for  a computer  required  consider 
able  care.  Our  declared  aim  was  to  reduce  to  negligible  proportions  the  numeri- 
cal errors  accompanying  the  implementation  so  that  consistent  results  can  be 
obtained.  This  has  been  achieved  at  very  moderate  extra  expense  in  terms  of 
computational  effort  and  thus  seems  well  worthwhile. 

The  programs  are  very  easy  to  use  and  have  been  applied  at  RAE,  Bedford,  to 
realistic  step  response  profiles  with  satisfactory  results  although  double 
precision  working  was  needed  because  of  the  shorter  wordlength  of  Che  SIGMA 
computer. 
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Scaling  and  shifting  the  response  to  a unit  step  gust 


Fig  3a&b  Numerical  approximation  of  a (smooth)  ramp  gust.  Tha  curvat  1.2.3 
ara  than  added  together  to  got  0(H.t).  These  diagrams  are  schematic;  the 
actual  computation  is  considerably  more  sophisticated  (see  text) 


ro  <i 


Fig  5 The  danger  of  sampling  the  step  response  (and  hence  of  generating  the  ramp 
response)  at  too  few  points.  The  primary  peak  and  overswing  will  both  be 
{ missed  because  the  programs  search  only  the  neighbourhood  of  the  extreme 

tabulated  points 


Fig  6 Variation  of  extreme  response  with  length  of  gust.  Defining  to  be  the 

envelope  of  the  above  curves  could  cause  trouble,^  the  worst  gust  pair 
would  be  wrongly  predicted  because  the  peak  at  would  be  ignored. 
Alternatively,  bemuse  y becomes  double  humped,  H might  be  missed 
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Fig  7a&b 


Step  response  having  the  properties 
li)  to3-to2>toi 

(ii)  area  A3  > area  A1  (A3b«A1) 


vtoi 

Vlt03-t02»  H 

b tytH)  for  the 

above 

step  response  using  the 

straight  ramp 

model 

Fig  7a&b  Showing  possible  discontinuous  behaviour  of  t^(H) 

Quadratic  intarpolation  of  a ikaw  function.  Points  1,  2,  3 are  the  starting 
values.  Points  4,  5,  6 are  obtained  by  fitting  a parabola  through  (1,2,  3), 
(2, 4,  3),  (4,  6,  3)  respectivaly  always  retaining  two  points  bracketing  the 
best  (highest).  This  process  is  inefficient  compared  to  that  described  in 
section  6 
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I AlUCRArT  RESPONSE  TO  SMOOTH  RaHP  GUSTS  (A.G.  PURCELL  APRIL  1976) 

y 


analytic  step  response  ObFlNEU  BY  ALPHA  > S.UO  AND  XI  ■ 0.50 


i aircraft  velocity  ■ 100.  FT/Str 

RAMP  RESPONSE  CALCULATED  AT  U.^0  SEC  INTERVALS  UP  TO  T > 10.00  SECS 

gusts  of  LFNGTM  h < ^0.0  FT  AHE  EFFECTIVELY  STfPS 

! SUCH  H(I>  IN  THE  DATA  ARE  IGNORED 


il 

1 

i 


Fig  9 Title  page  of  output  from  the  smooth  ramp  rasponie  program 


Fig  10 


CUBIC  SPLINE  FIT  TO  STEP  RESPONSE 
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Fig  10  Spline  fit  to  F(t) . Output  only  if  IPR  - 1 or  3 
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Fig  1 1 Values  of  the  integral  in  equation  (10)  at  intervals  of  At  = 0.2  s. 
Output  if  IPR  = 2 or  3 
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Fig  12  (contd) 
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